/*******************************************************************************

                                PERTURBHKL

  Title:        PERTURBHKL - Perturb fcalc file.
  Author:       Hanna Szoke
  Date:         6/23/97
  Function:     Using as input a model fcalc file, this program generates 
		complex vectors with Gaussian random perturbations to both
		amplitudes and phases.
	
*******************************************************************************/

#include "util.h"
#include "cellparams.h"
#include "symmetry.h"
#include "dims.h"
#include "sp.h"

float	ran1(int *) ;		/* packaged in util.c */

static	int	pseed ;			/* starting point for random numbers */
static	real	pert_frac ;

static	void	 getRpert(COMPLEX *, char *) ; 

void
perturbhkl_main(int argc, char *argv[])
{
     int	readFcList(char *) ;		/* in hklread.c */
     int	read_HDFfcalc(char *) ;		/* in hdfutil.c	*/
     int	HDFfilelength(char *) ;
     int	HDFdatalength(char *, char *) ;
     void       alloc_fc(int) ;

     char	fc_filename[MAXSTRING] ;
     char	out_filename[MAXSTRING] ;

     COMPLEX	*Rpert ;
     int	s = 1 ;

	/***********************************************************
	Check for right ballpark, etc.
	***********************************************************/

     if (argc > optind+5) 	/* too many arguments */
	  ballpark(caller) ; 
     else if (argc > optind+2) {
          sprintf(fc_filename, "%s", argv[optind+2]) ;
          pert_frac = atof(argv[optind+3]) ;
          if (argc == optind+5)
	       s = atoi(argv[optind+4]) ;
     } 
     else {
	  prompt("\nPlease enter filename, pert_frac, & optionally, seed: ") ;
          sscanf(terminp, "%s %f %s", fc_filename, &t1, message) ;
	  pert_frac = t1 ;
	  if ((int)strlen(message) > 0)
	       s = atoi(message) ;
	  sprintf(message, 
		"Running %s on file %s with pert_frac %f and seed %i", 
	         caller, fc_filename, pert_frac, s) ;
	  printTwice(message) ;
     }

     hello(caller) ;

     if (s < 1) 
	 Fatal("Bad seed - stopping.") ;
					    
	/***********************************************************
        Just to ensure that we don't use very small seeds that have 
	the weird property that they occasionally return 1, add a 
	"random" number before complementing.
	***********************************************************/
							    
     pseed = -(s + 137) ;
							  
     sprintf(out_filename, "pert%02d.fcalc", s) ;

     readBasicInput() ;

	/***********************************************************
        Allocations and reads.
	***********************************************************/

     if (HDFIN)
          Nfc = HDFdatalength(fc_filename, "fcalc/Nfc") ;
     else
          Nfc = filelength(fc_filename) ;

     alloc_fc(Nfc) ;

     if (HDFIN)
          Nfc = read_HDFfcalc(fc_filename) ;
     else
          Nfc = readFcList(fc_filename) ;

     Rpert = (COMPLEX *) e_malloc(Nfc*sizeof(COMPLEX), "perturbhkl_main") ;
     getRpert(Rpert, out_filename) ;
}

				/****************************************
				****************************************/

void	 getRpert(COMPLEX *Rpert, char *filename) 
{
     real	normGauss(int *) ;
     real	cleanPhase(real *) ;
     FILE	*fpout ;

     COMPLEX	R ;
     real	omega ;
     real	amp0 ;		
     real	Ramp, Rphase ;
     real	newamp ;
     int	n ;

     fpout = fopen(filename, "w") ;
     if (fpout == NULL)
     {
          sprintf(message, "Cannot open %s", filename);
          Fatal(message) ;
     }

     amp0 = *(Fc_relist) ;
sprintf(message, "amp0 = %g\n", amp0) ; 
printTwice(message) ;

     for (n = 0; n < Nfc; n++) {

          /*** next 12 lines are equivalent of getAmpPhase() ***/
          R.re = *(Fc_relist + n) ;
          R.im = *(Fc_imlist + n) ;
          Ramp = sqrt(R.re*R.re + R.im*R.im) ; 

          if (Ramp < amp0*AMP_EPSI)  {
               Ramp = 0 ;
               Rphase = 0 ;
          }
          else {
               omega = RTOD * atan2(R.im, R.re) ;
               Rphase = cleanPhase(&omega) ;
          }

	  newamp = Ramp + normGauss(&pseed)*pert_frac*Ramp ;

          fprintf(fpout, "INDEX %10.4f %10.4f %10.4f FCALC   %g   %g\n",
               *(Fc_Hlist+n), *(Fc_Klist+n), *(Fc_Llist+n), newamp, Rphase) ;
     }
     fclose(fpout) ;
}
